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Abstract 

> 

The production of acoustic signals from the interactions of ultra-high energy (UHE) 
cosmic ray neutrinos in water and ice has been studied. A new computationally fast 
and efficient method of deriving the signal is presented. This method allows the 
implementation of up to date parameterisations of acoustic attenuation in sea water 
and ice that now includes the effects of complex attenuation, where appropriate. 
On The methods presented here have been used to compute and study the properties 

of the acoustic signals which would be expected from such interactions. A matrix 
method of parameterising the signals, which includes the expected fluctuations, is 
also presented. These methods are used to generate the expected signals that would 
be detected in acoustic UHE neutrino telescopes. 
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1 Introduction 



In recent years interest has grown in the detection of very high energy cosmic 
ray neutrinos which offer an unexplored window on the Universe pQ. Such 
particles may be produced in the cosmic particle accelerators which make 
the charged primaries or they could be produced by the interactions of the 
primaries with the Cosmic Microwave Background, the so called GZK effect 
[2] . The flux of neutrinos expected from these two sources has been calculated 
[3f4] . This is found to be very low so that large targets are needed for a 
measurable detection rate. It is interesting to measure this neutrino flux to 
see if it is compatible with the values expected from these sources, with any 
incompatibility implying new physics. 

Searches for cosmic ray neutrinos are ongoing in AMANDA pj], IceCube [6], 
ANTARES [7], NESTOR 0, NEMO [9], KM3NeT pE0] and at Lake Baikal 
[TT] detecting upward going muons from the Cherenkov light in either ice or 
water. In general, these experiments are sensitive to lower energies than dis- 
cussed here since the Earth becomes opaque to neutrinos at very high energies. 
The experiments could detect almost horizontal higher energy neutrinos but 
have limited target volume due to the attenuation of the light signal in the 
media. The Pierre Auger collaboration, using an extended air shower array de- 
tector, are searching for upward and almost horizontal showers from neutrino 
interactions [12J. In addition to these detectors there are ongoing experiments 
to detect the neutrino interactions by either radio or acoustic emissions from 
the resulting particle showers pQ. These latter techniques, with much longer 
attenuation lengths, allow very large target volumes utilising either large ice 
fields or dry salt domes for radio or ice fields and the oceans for the acoustic 
technique. 

In order to test the feasibility of detecting such neutrinos by the acoustic tech- 
nique it is necessary to understand the production, propagation and detection 
of the acoustic signal from the shower induced by an interacting neutrino in a 
medium. This has been treated in some detail in [T3j, however, in this treat- 
ment it is difficult to incorporate the true attenuation of the sound which has 
been found to be complex in nature [14] in media such water. Such 

complex attenuation causes dispersion of the acoustic signal and complicates 
both the propagation of the sound through the water and the signal shape at 
the detectors. 

This paper is organised as follows. In section 2 the new approach to calcu- 
lating the acoustic signal pressure is described and section 3 describes the 
methods used to model the attenuation of the sound as it propagates through 
the medium. Section 4 describes the detailed calculations of the sound sig- 
nal in water and in ice as it arrives at the detector starting from the shower 
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simulations described in [15]. Finally a new method of simulating signals in- 
corporating shower-to-shower fluctuations is described in Section 5. 



2 The thermo-acoustic mechanism 

The standard equations used to determine the thermo-acoustic integrals are 
outlined in [13]. In this paper we use a complementary approach. 

For the thermo-acoustic mechanism even though it is the pressure, p, that is 
detected it is the volume change, Q s , which couples to the velocity potential, 
$, which in turn creates the sound wave. P] Interestingly, the velocity potential 
as a concept precedes the magnetic vector potential by over 100 years and was 
introduced by Euler in 1752 [PT] . 

Three of the most important variables in acoustics are the pressure change 
from equilibrium, p, the particle velocity, v, and the velocity potential, $. 
Assuming zero curl these three variables are related by: 



where p is the density. For sources of acoustic energy this velocity potential 
has a function in acoustics equivalent to the magnetic vector potential in elec- 
tromagnetism. We are trying to solced the wave equation to get the pressure 
pulse at the location of an observer placed at \r\. This can be dome by in- 
tegrating all the contributions of infinitessimally small sources at locations 



For an observer at r and a shower event at r' separated by a distance r = 



where q s is the time rate of change of an infinitesimal volume and c the velocity 
of sound. In our case we are interested in integrating q s over the cascade volume 
where this volume change is caused by the injection of an energy density E (in 



2 This is analogous to the situation in the radio detection of UHE neutrinos where 
even though it is the electric field E that is detected, it is the current density, 
j, which couples to the magnetic vector potential, A, which in turn creates the 
electromagnetic wave [16]. 



v = V$ 



p = -p 



<9$ 



(1) 




(2) 
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J.m 3 .s 1 ) over this cascade volume. For an infinitesimal volume the volume 
change starting at time t = 0, is given by: 

q s {¥, t) = 4rE(F, t) - 1 / q s (7, f)dt', (3) 



pC p T 







where (3 the thermal expansion coefficient, C p the specific heat capacity, p the 
density and r the thermal time constant. The integral term is caused by cooling 
as the deposited energy, within the volume, conducts or convects away into 
the surrounding fluid. However as the time constant for this thermal cooling 
mechanism is of the order of tens of milliseconds [19] , and as we are primarily 
interested in the case where the energy is injected nearly instantaneously in 
acoustic terms (ns), the second term in equation 3 can be ignored as it is 
about six orders of magnitude lower than the first term. Equation 3 can then 
be integrated over the entire volume cascade using Green's functions (see for 
example [18]) and yields for an observer at distance, r, from the source: 

= JJ ^^M¥,t') {t=t , + ^ /c) dVdt' (4) 
v ' ' 



where t — r/c is the retarded time, this simply implies that the contribution 
from each point in the source, in both space and time, travels to the listener 
with the speed of sound. Equation 4 lends itself to efficient numerical solution. 
If the energy deposition is modelled using Monte Carlo points with a density 
proportional to energy, will be a scaled histogram of the flight times to 
the observer. The pressure can now be derived from equation 1. A further 
simplification can be made if the energy deposition as a function of time 
is identical at all points in the volume. The velocity potential can then be 
calculated from the convolution integral: 

^'ikpl ' j v \KFWE t (i-X)^ (5) 



where A is a dummy variable used to evaluate the integral. This can be sim- 
plified further since E(t) can be approximated by 5(t) for the case we are 
interested in. The time dependence on the right hand side will disappear as 
the integral of a delta function is equal to one. This formulation of the thermo- 
acoustic mechanism leads to a solution in the far field which is proportional 
to the Griineisen coefficient 7g |2Q] . If, after an arbitrary 3D rotation, the 
observer is, for example, at a distance from the source along the x' axis much 
further away than the dimensions of the source, then with increasing observer 
distance, x' , the E(y') and E(z') contribution will more and more closely 
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resemble delta functions as the spread in arrival times caused by their contri- 
butions approaches zero. Equation 5 reduces to: 



d>(t) 



AirCpp x' — x' 



7 E(c(t-t )) 



(3 c 
AnC p p x'q 



E(c(t-t )), 



(6) 



where to is the flight time from the centre of the shower to the observer. 
The velocity potential, $(t), is simply a scaled cross section of the energy 
deposition and pit) a projection of $(t): 



The (3c 2 /C p term is the Griineisen coefficient and gives the relative acoustic 
pulse heights for different media. With this formulation the far field solution 
in the absence of attenuation makes a number of predictions: 

a) The pressure pulse is a scaled derivative of the projection along the line 
of sight to the observer; 

b) Only the distribution along the line of sight is important. Hence, for 
example, with an observer on the 1 J deposition into a Gaussian 
sphere (i.e. where x,y and z are randomly distributed with a normal 
distribution with the same standard deviation, (J X ^ )Z ) with a = 1 cm, 
centred on the origin, will observe the same pulse as a tri-axial Gaussian 
distribution (as above where the standard deviations in x, y and z are 
different) with o x = 1 cm, a y = 1 mm and o z = 1 m centred on the 
origin. This is the case when x » <J y ,a z ; 

c) The amplitude of the pressure pulse will depend on 1/cr 2 where o"o is 
the projected standard deviation in the acoustic pulse height as seen 
from the observer. Hence an observer on the y axis will see a pulse one 
hundred times greater in magnitude than one on the x axis for the tri- 
axial distribution as described in item b) above; 

d) At angles greater than a few degrees from the x — y plane the amplitude 
of the pressure pulse will depend on 1/ sin 2 9. 

This is illustrated in figure [T] where each of the sub- figures illustrate the cor- 
responding point in the list above. 




f t E(c(t-t )) 



(7) 
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3 Attenuation 



3. 1 Introduction 



The acoustic pulse is affected by the medium through which it travels and it 
acts as a filter causing frequency dependent attenuation, a(u). Attenuation of 
acoustic pulses is caused by a combination of absorption and back scattering. 
The amplitude is attenuated by a factor given by: 

a(u) = e- fc(w)r , (8) 



where r is the distance travelled and k is the attenuation coefficient (Nepers/m). 
The resultant pulse can be determined by converting the pulse into the fre- 
quency domain by taking a Fourier transform, multiplying and taking the 
inverse Fourier transform: 

oo / oo \ 

*«(0 = ^ / «H / Ht')e-^dt) e^dw, (9) 

— oo Voo / 



where <£> a is the velocity potential with attenuation and t' = t — r/c. Alter- 
natively the effect of attenuation can be determined by taking the inverse 
Fourier transform of the frequency characteristic (i.e. the modification of the 
pulse between transmission and reception due to the medium) and convolving 
this with the un-attenuated pulse: 

oo 

a(t) = I J a{u)e^du 

-OO /1Q\ 

OO oo \ u / 

$ a (t')= / a (t' - \)®(X)dw = J $(t! - X)a(X)d\, 

— oo — oo 



where A is a dummy variable used to evaluate the integral and a(t) is the pulse 
created by passing a Dirac delta function source through the medium and is 
referred to as the unit impulse response The equivalence of these methods is a 
statement of the convolution theorem: multiplication in the frequency domain 
is convolution in the time domain. Values of the pressure attenuation length in 
sea water (ice) are 3.8 km (10 km) at 10 kHz and 2.6 km (8.5 km) at 25 kHz. 
Attenuation lengths for intensity are a factor of two less. 
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3.2 Attenuation in distilled water 



Evaluating a(t) analytically is normally difficult. However for a liquid in which 
the attenuation is dominated by viscosity it is straight forward. This is the 
case for example in distilled water. It can be shown (see, e.g. [13J) that in 
this case the attenuation coefficient, k, and attenuation, a, for frequencies of 
interest, are given by: 

k = — — , aiuf) = e _2 = V = e 2 <m^) 2 where aJuo) = \[^-, (11) 
2c uq V r 



and uiq is medium dependent. The attenuation a(uj) is Gaussian in shape with 
a standard deviation of a a which is proportional to In distilled water, 

ujq ~ 10 12 rad s^ 1 and Lehtinen et al. [2T] have used a value of 2.5 x 10 10 rad 
s _1 to approximate the attenuation of tropical sea water in the 10 4 — 10 5 Hz 
region. The unit impulse response is given by the inverse Fourier Transform 
as: 



a(t) = i / e-^ 2 e^du = J Q^e 



-oo 

,2 



_ (12) 

e 2 "a(t) 2 where cx a (t) 



and is also a Gaussian distribution with a standard deviation which gets wider 
with the square root of the distance. Indeed a a {t) = l/a a {uj). If the source term 
is also Gaussian the overall pulse profile can be evaluated. The convolution 
integral of two Gaussian distributions yields a Gaussian distribution with the 
standard deviations adding in quadrature. Therefore if: 



f3E 1_ 



•I- 



) = where a, a = Jal + al 



(13) 



then the pressure p is given by: 



This is equivalent to equation 18 in reference [T3], but is derived using a quite 
different approach. 
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3.3 Attenuation in sea water and ice 



Acoustic attenuation in seawater is almost totally caused by absorption. In 
the l-100kHz region it is dominated by a chemical relaxation process that is 
connected with the association-dissociation of magnesium sulphate (MgSO^) 
ions under the pressure of the sound wave. Below 1kHz a similar mechanism 
involving boric acid (B(OH) 3 ) is responsible for much of the observed atten- 
uation. Taken together these mechanisms result in an attenuation of acoustic 
waves in seawater and a velocity of sound which are both frequency depen- 
dent. Experimentally however whereas it is straightforward to estimate the 
magnitude of \z\ of the complex attenuation it is very difficult to determine 
the phase angle Zz and no current measurements exist in the literature. 

The magnitude of the attenuation however is well measured and the definitive 
work in this area is considered to be that of Francois and Garrison [22] . More 
recently, Ainslie and McColm |23j have published a simplified parameterisa- 
tion of the magnitude of the attenuation in sea water which maintains a similar 
accuracy to the parameterisations of Francois and Garrison. This is a function 
not only of frequency but also depends on depth, z, salinity, S, temperature, T 
and pH. Whereas there are no direct measurements of the attenuation angle, 
Lieberman [21] gives a clear presentation of the chemical processes causing 
the attenuation while Niess and Bertin [25] have published a complex attenu- 
ation formula based on Mediterranean conditions. Here we present a complex 
version of the Ainslie and McColm formulation, which retains the attenuation 
magnitude but introduces the phase shifts predicted by Lieberman. In essence 
the attenuation consists of three components, two of these are complex, high 
pass filters with cut off frequencies ujb (~27r x 10 3 rad s -1 ) for boric acid and 
u>M g (~27r x 10 5 rad s _1 ) for magnesium sulphate. The third is the pure water 
component, which is real. The ACORNE parameterisation uses a x values that 
are the respective attenuation coefficients in dB/km: 



lu b = 15607iy§ e T / 26 , uo Mg = 84000vre T / 17 

a B ~ 2n e 56 , a Mg ~ ^ (1 + ^) 



(15) 



_ 49xl0~ 9 n -(T/27+z/17) 



adB/km - + s+LUMg + a w u wnere s-iu. 



In ice the mechanisms are less well understood. However following Price 
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in regions of deep ice of most interest for acoustic neutrino detection, e.g. the 
South Pole, the dominant attenuation mechanisms below a few hundred kHz 
are: absorption due to proton reorientation (relaxation) and scattering due to 
bubbles and grain boundaries. For South Pole conditions Price predicts the 
absorption length (apart from very low frequencies <100Hz) to be constant: 

r 2 

a L [m = — 16 



where 5 max is the log decrement at asymptotically high frequencies, r m is a 
temperature dependent relaxation time and is the velocity of sound in 
ice (3920 m s _1 ). The parameters are normally evaluated experimentally. For 
South Pole conditions these predict an absorption length of ~10 ±3 km at 
frequencies above 100 Hz. The dominant scattering mechanism is caused by 
grain boundaries and for the grain size expected at the South Pole the ice will 
behave as a Rayleigh medium. Scattering is proportional to the 4th power of 
frequency, up to ~300 kHz: 



6 i (* m -')=5xl0*f^-V(^y 
v ; V0.2cm J \ 10kHz J 



(17) 



where d is the mean grain diameter which is 0.2 cm. As the acoustic pulses from 
neutrino interactions tend to be highly directional we make the assumption 
that the attenuation is given by the sum of the absorption plus total scattering 
rather than the more usual assumption that the attenuation is given by the 
sum of absorption plus back-scattering. 

A comparison of the ACORNE attenuation parameterisation and earlier work 
is presented in figure |2j together with the anticipated attenuation in ice. The 
f 2 approximation assumes the sea behaves like a more viscous form of distilled 
water by decreasing the value of u to match the measured attenuation in the 
10 4 to 10 5 Hz frequency region for tropical waters, consistent with the location 
of the SAUND array [21]. The other curves are optimised for Mediterranean 
conditions (T = 15°C; S = 37; pH = 7.9; z = 2km). Since the Ainslie & 
McColm/ACORNE and Francois & Garrison results match very closely they 
are depicted as the same curve in the figure. Phase shifts for a pulse 1 km 
distant from the neutrino interaction are shown for the complex attenuation 
case, corresponding to a velocity increase of ~0.005% between 1 and 100 kHz. 
For ice the attenuation is constant up to a frequency ~25 kHz where Rayleigh 
scattering starts to dominate. The absorption in ice is lower than that of sea 
water in the important 10-100 kHz region. 

The effect of the various forms of attenuation on a point source of energy 
5(r)5(t) observed at 1 km from the origin is shown in figure [3j The three real 
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attenuation mechanisms in water give similar symmetric pulses. This gives 
confidence that the magnitude of the attenuation need not be modelled in 
fine detail. Interestingly however whereas the magnitude of the two complex 
attenuation models is similar to the other water models, (and, indeed identical 
in the case of the ACORNE and Ainslie & McColm models) the phase shifts 
cause considerable pulse distortion from the anticipated symmetric shape. This 
is a result of the non-linear nature of the phase shifts creating the group 
delay, dtp/duj, to vary, hence creating dispersion. The phase shift also has the 
effect of marginally reducing the pulse velocity at lower frequencies hence 
slightly delaying the pulse. In ice the pulse is again symmetric and is larger in 
amplitude by a factor of three. This is not, as one might expect, due to the 
larger value of the Griineisen coefficient in ice, but is a direct consequence of 
the lower attenuation in the 10-100 kHz region. []j The ripples before and after 
the main pulse are caused by the abrupt f 4 scattering dependence switching 
on yielding a rapid rise in attenuation. Any sharp transition in the frequency 
domain will tend to cause ripples in the time domain, these ripples are a real 
physical effect and not a numerical artefact. 



4 Acoustic integrals for CORSIKA simulated pulses in water and 
ice 

4-1 Introduction 

In this section we compute full 3D simulations of acoustic pulses from our pre- 
viously published neutrino shower simulations [IS]. In our previous work, 1500 
showers were modelled using the air shower simulation program CORSIKA [27] 
modified for a water or ice medium [15] corresponding to 100 showers in each 
half decade increment from 10 5 to 10 12 GeV. The deposited energy was binned 
in a cylinder of 20 m in length and 1.1 m in radius. Longitudinally 100 bins of 
width 20 cm with bin centres from 10 cm to 1990 cm were used. Radially, 20 
bins were used: ten of width 1 cm with centres from 0.5 to 9.5 cm and 10 of 
width 10 cm with centres from 15 cm to 105 cm. The contribution of the outer 
ten bins towards the acoustic pulse is minimal; they are included to capture all 
the energy in the shower and are useful for conservation of energy checks. This 
spacing should be sufficient to accurately model the acoustic frequency com- 
ponents in the pulse on axis up to ~75 kHz in water and ~200 kHz in ice. In 
this section the average shower distribution for each of the 15 half decades in 
energy is used directly as an input to drive the acoustic integral. The neutrino 

3 The ratio of (3/C p is 1.38 between ice and water and the 1/c 2 term is a factor of 
6.83 larger for ice. However, the 1/c 2 term only applies to extended sources and can 
be ignored for a point source. 
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creates a hadronic shower of approximately 10 m in length and 5 cm in radius. 
The thermal energy is effectively deposited instantaneously in acoustic terms, 
creating an acoustic radiator analogous to a broad slit. In the far field, using 
Fraunhofer diffraction, the angular spread of the radiation is approximately 
X/d. The wavelength A is of the order of the diameter of the cylinder and d 
the length of the cylinder, yielding a pancake width in the order of 1 degree. 
Water and ice are different media, the suitability of the medium for acoustic 
detection will depend largely on the magnitude of the created acoustic pulse 
and the attenuation. For extended sources of this nature a suitable figure of 
merit to compare different media is given by the ratio of the product of the 
Griineisen coefficient (7c) and the anticipated energy deposition (dE/dX): 



Using this expression, ice should produce acoustic pulses which are about an 
order of magnitude greater in amplitude than water. However because the 
velocity of sound in ice is about 2.5 times that of water, the frequencies in 
the pulse will also be about 2.5 times higher causing ice to have a greater 
attenuation than water. 



4-2 Coordinate system 

In the subsequent analysis a coordinate system is chosen such that the neutrino 
interacts at the origin and travels vertically along the z axis where the value of 
z increases with depth and the origin. The point (0, 0, — z p ) is chosen such that 
the maximum "pancake" energy (i.e. acoustic, and not deposited energy) at 
1 km from the shower is at zero degrees. The value of z p is energy dependent 
varying from approximately 460 cm at 10 5 GeV to 780 cm at 10 12 GeV. If the 
radial cross-section of the shower were constant one would expect the pancake 
maximum to correspond to the centroid of the energy deposition taken along 
the longitudinal axis. In reality however the radial distribution gets broader 
with the age (depth) of the shower; the earlier part contributes more to the 
pulse energy than the later part of the shower. The point (0,0,-z p ) can be 
determined in two ways, either from the position of the shower centroid, z c or 
from E, the primary energy (in GeV) which defines the pancake plane. A fit 
to sea water data yields following relationships for z p : 



e = 



1GP 



dE c 2 (3 dE 



(18) 



z p = 1.05^-0.874 

z p = -0.046 (log 10 E) 2 + 1.3 log 



E - 0.84. 



(19) 
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In Antarctic ice the pancake depth has to be increased by approximately 6% 
due to the relative densities of the media. 

4-3 Acoustic integrals in sea water 

For the initial analysis the observer is positioned at 1 km from the shower in 
the centre of the acoustic pancake to allow easy comparison with previously 
published results [21] [25]. Complex attenuation (equations 15) was assumed 
and the acoustic integral calculated for each of the 15 half decades in energy 
as described above. The results are plotted in figure [4} As the integrals have 
been calculated with high precision, error-bars are not plotted, as they are less 
than a line width. As can be seen a characteristic bipolar pulse is produced. 
In figure |4^i) the pulse shape is plotted for three energies. The maximum pulse 
amplitude is normalised by energy. As can be seen the pulse shape is very 
similar for the three chosen energies (10 5 , 10 8 ' 5 and 10 12 GeV). The pulse 
height however seems to scale slightly more rapidly than energy. In figure |4Jd) 
the maximum and minimum pulse amplitudes are plotted on a log-log plot. 
The fitted lines are constrained to be proportional to energy. It is clear from 
figure |4]d) that any increase in pressure over proportionality is minimal. A 
linear fit yields: 

log 10 (Ana X ) = 1.0021(±0.0018) log 10 (£) - 11.93(±0.0012), (20) 

where p ma x (Pa) is the maximum pressure and E the energy in GeV. The 
errors quoted are the statistical errors from the fit. This yields the result that, 
to a good approximation, the maximum pulse height at 1 km is 1.22 pPa per 
GeV in the plane of the pancake. The average frequency (using the MUSIC 
algorithm [28]) is very stable with energy and is approximately 26 kHz. The 
asymmetry (\P max \ - \P m in\) / {\Pmax\ + \Pmm\) is similarly stable and is about 
0.2. 

Consider now figure [5] showing the angular spread of the pancake. As the 
acoustic pulse will most likely be detected by a matched filter (see, for example, 
[21] for a full discussion of matched filters), which integrates over the pulse, 
it is the integrated pressure, or square root of the pulse energy which is of 
most interest. This is plotted in figure [5^,), which shows the angular spread 
of the shower from —4° to 4°. The pancake is extremely narrow and narrows 
further with increasing energy. In figure |5Jd) the full width of the pancake at a 
pressure levels of 50% and 10% of maximum is illustrated. As anticipated from 
simple far field diffraction theory (the diffraction minimum is at 8=\/d) the 
spread of the pancake decreases with increasing shower energy. This is largely 
because the shower gets longer as the energy increases causing a narrowing of 
the pancake. 
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It is interesting to cross check the acoustic mechanism by looking at the energy 
flowing through a 1 km integrating sphere. The energy of the acoustic pulse 
flowing through each square metre of surface (the fluence) is given by: 



CO 



P 2 (t,r,6) 
Z 




dt, 



(21) 



— oo 



where Zq is the characteristic impedance (~1.5 x 10 6 kg m~ 2 s" 1 for water 
and ~3.7 x 10 6 kg m -2 s _1 for ice) and 9 is the angle out of the plane of the 
pancake, the radiation is cylindrically symmetrical. This can be integrated 
over a sphere, in figure [5b) the result of this integral is displayed. The fit of 



captured acoustic energy, E c vs. deposited energy E (see figure |5ji) )is given 




where k= 22.8 without attenuation and 23.1 with attenuation. Hence at 1 km 
attenuation reduces the acoustic energy to 50% of its value without attenu- 
ation. The linear coefficient is 2.00 to within the statistical accuracy of the 
showers generated by CORSIKA, indicating that the efficiency rises quadrat- 
ically with energy. Due to the coherent nature of the acoustic emission mech- 
anism, the acoustic pulse amplitude depends linearly and the acoustic energy 
depends on the square of the deposited energy (assuming constant shower 
shape). This coherent behaviour breaks down at energies far beyond those of 
interest here. 

It is interesting to look at the effect of complex attenuation on the pulse asym- 
metry in water as there is no guarantee that the pulse will become symmetrical 
even in very far field. Figure [6] illustrates the mean pulse frequency, pressure 
times distance and asymmetry for a 10 11 GeV energy deposition and for an 
observer in the pancake plane for distances of 10 m to 10 3 km. In the absence 
of attenuation, once far field has been established (~1 km) both p max d and the 
mean frequency, f max should be constant and the asymmetry zero. In practice 
the mean frequency drops with distance, as higher frequencies are more quickly 
absorbed than lower frequencies, falling from ~50 kHz close to the shower to 
~200 Hz at 10 3 km. The product of maximum pressure and distance rises ini- 
tially as the radiator approaches far field conditions. The maximum pressure 
then drops both because the energy is absorbed and the pulse becomes more 
spread out in time. The asymmetry starts with a value of 0.7 dominated by 
near field effects then drops to ~0.1 at 1-10 km and has a peak of ~0.5 at 100 
km. The dispersion will be a maximum at frequencies around the resonant 
peaks of MgS0 4 and B(OH) 3 . The B(OH) 3 peak is at ~1.3 kHz causing the 
rise in asymmetry where the mean frequency matches this value. 




by: 



log 10 £ c = 2.001og 10 £ 



-k, 



(22) 
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In figure [7] the mean frequency, maximum and minimum pressures and asym- 
metry are plotted versus angle and distance for a 10 11 GeV shower using the 
ACORNE complex attenuation (T=15°C, S=37ppt, pH=7.9, z=2km). As the 
neutrino is travelling vertically downwards then positive angles (left hand side 
in figure [7]) are measured out of the plane of the pancake and towards the sur- 
face of the water. In figures [7^,) and [7Jd) the maximum and minimum values 
of P in dB re. 1 Pa are displayed. As anticipated, the pressure decreases with 
angle and frequency. Figure ^f)p) shows the asymmetry where two effects are 
particularly noteworthy. Firstly the dominant effect of complex attenuation 
drives the pulse towards positive asymmetry. Secondly there are two regions 
with asymmetry of greater than 0.6. In these regions the geometry causes a 
spike in the velocity potential yielding a very non bipolar pulse. Consider now 
figure showing the mean frequency. As can be seen at angles below ~1° 
and distances below 100 m the frequency with maximum energy is above 40 
kHz. This reduces with both angle and distance. The decrease of frequency 
with angle is caused mainly by geometric projection as discussed in section 2. 
The decrease of frequency with distance is caused mainly by high frequencies 
being attenuated more rapidly than low frequencies. 



4-4 Acoustic integrals in ice 



Following the procedure outlined in section 3.2 and adopting a model of 
Antarctic ice from [26] , the initial analysis was again at 1 km from the cascade 
in the plane of the pancake. (The pancake is at a depth about 5% further from 
the shower origin due to the relatively lower density of Antarctic ice to that 
of sea water). In figure [8^,) the detected acoustic pulse is plotted for three 
energies. The pulse shape again scales with energy but is 5-6 pPa per GeV, 
(about 5 times that of water). The pulse is narrower than in water, though not 
by the factor of 2.6 as predicted by the relative velocities. The pulse is more 
symmetric as the attenuation is dominated almost entirely by scattering and is 
non-complex. The pulse also shows a ripple indicative of the sharp f 4 nature of 
Rayleigh scattering. In figure^) the relationship between the maximum pulse 
height and energy is displayed. Again the pulse height grows almost linearly 
with energy and the best straight line fit of log 10 p to log 10 E as in equation 
20 yields a slope of 1.0079 (± 0.0017) and an intercept of -11.349 (± 0.151), 
corresponding again to a growth in pulse height which is very slightly greater 
than energy deposition and an intercept corresponding to 4.48 pPa/GeV. In 
practice it is the higher energies which are of most interest. If the fit is done 
in the 10 9 -10 12 GeV region the intercept corresponds to 5.35 pPa/GeV. The 
average frequency is nearly constant at 39 kHz and the asymmetry is nearly 
zero. 

Consider now figure [9] showing the angular spread of the pancake in ice. The 
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analysis is identical to that described in section 4.3. As can be seen in fig- 
ures |9^i) and^), the angular spread of the pulse is very similar to that of sea 
water, but is slightly broader due to the effective wavelength in ice being longer 
than that of water. In figures [9J:) and |9]l) the energy captured on an integrat- 
ing sphere at 1 km is calculated. In figure [9]i) the scaling of captured energy 
with the square of the deposited energy is again evident. The intercepts are 
-22.19 GeV and -21.76 GeV with and without attenuation respectively. Hence 
about 11 times as much energy is created in ice as in water. Once attenuation 
is included the energy drops to 37%; 63% of the energy is lost nearly entirely 
due to scattering. This loss is slightly more than water largely because ice 
pulses are at a higher frequency causing more effective attenuation. 

As with sea water the analysis is extended to distances between 10 m and 
50 km and angles between ±90°. The mean frequency, maximum and mini- 



mum pulse heights and asymmetry are displayed in figure 10 Here the mean 
frequency is defined as the frequency of a single cycle sine wave that most 
closely represents the acoustic bipolar pulse, this is usually a little higher 
than the peak frequency. As anticipated both the frequencies and pulse height 
are greater than that of sea water. Also, due to the real nature of the atten- 
uation the asymmetry falls rapidly to zero in the plane of the pancake and 
the odd symmetry in the far field between plus and minus angles is also more 
evident. 



5 Modelling showers and parameterising the fluctuations 



5. 1 Introduction 



Previous work in this area has [15] (25] [21] concentrated on modelling average 
shower parameters as a function of energy. Fluctuations have not been con- 
sidered. Indeed due to the nature of previous parameterisations which involve 
correlated variables, the inclusion of fluctuations is not feasible as varying one 
parameter will mean that the other parameters have also to be tuned. The 
strategy in this section is to parameterise the shower and its fluctuations using 
an orthogonal basis set. 



5.2 CORSIKA Monte Carlo data 



The showers have been modelled with CORSIKA [27J, which uses a thinning 
process. The stochastic fluctuations in individual showers were smoothed us- 
ing a non-causal 3rd order Butterworth filter. In figure 11 the effect of the 
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Butterworth filtering is shown for four typical 10 5 GeV showers. The overall 
shape of the shower is retained but the noise considerably reduced. 



5.3 SVD parameterisation 



Parameterisation using singular value decomposition (SVD) |29j is an eigen- 
vector based technique and has a number of advantages which include appli- 
cability to data for which functional parameterisation is difficult and, most 
importantly, the ability to include fluctuations. The ability to include fluc- 
tuations stems from the orthogonal nature of the parameterisation; each pa- 
rameter can be varied independently as there is no covariance. The standard 
SVD method is matrix based and directly applicable in a 2d parameterisation. 
In the case of a 3d parameterisation the SVD has to be applied a number of 
times; we use two in this study. As the radial distribution is most critical 
for the acoustic pulse shape, this is used as our primary dimension. Each of 
our 1500 CORSIKA generated showers Sk is treated as a matrix with 100 
rows corresponding to longitudinal distance and 20 columns corresponding to 
radius. These are appended to form: 



O 



V Sk 7 



(23) 



creating an observation matrix O of size 150,000x20. A singular value decom- 
position was performed on the O matrix: 



O = WLV = (W S W N ) I 8 ° I I " " I , (24) 



L 



N 




where W and V are unitary containing the column and row eigenvectors of 
OO t respectively sorted in decreasing magnitude of the associated eigenvalue 
and L is diagonal containing the square roots of the respective eigenvalues. 
The S and N subscripts correspond to signal and noise and are discussed 
below. If all 20 eigenvectors are used the data can be reproduced perfectly. 
With fewer eigenvectors the data can be approximated with an accuracy a 
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given by: 



a = 



J2 Lj, 

1 = 1 

i=20 
i=l 



(25) 



The matrix has been partitioned into a signal space and noise space by choos- 
ing the n eigenvectors which contain the most information and assigning those 
to the signal space, the remaining 20 — n eigenvectors become the noise space: 



O = ()s On W s I LS Vs+W N ° 

J I L N 



V n . 



(26) 



The choice of the value of n can be somewhat subjective, however choosing 
n = 4 yields a = 0.95, which is more than sufficient given the other errors in 
the process. The matrix V s (20x4) consists of the four row eigenvectors. The 
matrix W s has four columns (4x150,000), this can be partitioned as follows: 

W 21 W 3 i W 41 ^ 



w 12 w 22 w 32 w 42 



V 



w 1k w 2W w 3k w 



4k 



(27) 



where each of the sub matrices W mn is of size 100x1. However these have 
the shower information sequentially and need to be reshaped to look at the 
longitudinal correlations between showers. The matrix is reordered into four 
matrices of size 1500x100: 



W s = 



'w20 
w? 2 



v w? k/ 



wj 2 



(28) 



\W 4 T k y 



A further singular value decomposition is applied to each of these matrices in 
turn. As these matrices contain successively less information four of the column 
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eigenvectors were retained for Wi, three for W2, two for W3 and one for W4. 
The four L m and V m matrices of this second SVD contain the information 
yielding the 10 parameter values for each shower. These parameters may be 
recovered by creating a matrix: 



V'i V' 2 V' 3 V' 4 



( L'i 








\ 





L 2 














L's 





V 








L' 4y 



(29) 



The matrix A (1500x10) contains on each row the coefficients a\ to aio for 
each shower in turn. The original showers, Sk, can simply be reconstructed by 
matrix multiplication of the appropriate coefficients: 

S k = a kA W 11 V 1 + •• + a(c,4Wi 4 Vi + afc,5W 2 iV 2 + . 
+afc,7W 23 V2 + a fe)8 W 3 iV 3 + a fci9 W 32 V 3 + a M0 W 4 iV 4 . 



The showers are recreated to within an average accuracy of 5 % in profile 
and magnitude. Lower energy showers have the greatest variation in shape. In 
figure [TT] the longitudinal distribution of four of the most extreme showers are 
plotted. These are chosen by having the greatest and least component of the 
first two radial eigenvectors. These will be some of the most difficult showers 
to model. The reconstructed showers though clearly not perfect do however 
reasonably fit the data. The mean and standard deviation of the parameters 
in each of the 15 half decades can now be used to reproduce the statistics 
of the showers, a\ to a 10 . Interpolation can be used to model intermediate 
values. Care must be taken however as these are not orthogonal as is normally 
the case with the SVD. Because of the linear nature of the parameterisation 
however these correlations can be reintroduced by multiplying by the matrix 
square root of the correlation matrix: 

A(E) = A(E) + ^RiA^E)^ C(E) 1 *, (31) 



where A(E) is a vector containing the average aj values at a given energy, R 
is a vector of normally distributed random numbers with a mean of and a 
standard deviation of 1, A a (E) are the standard deviations of ai and C(E) is 
the matrix of the correlation coefficients between parameters within a shower 
and averaged over all showers at a particular energy. The parameter values 
are too cumbersome to list here but are available on the ACORNE web site 
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5.4 Fluctuations applied to acoustic pulses 



In figure 12 the anticipated pulses from 100 CORSIKA generated events for 
each of the energies 10 5 , 10 7 ' 5 , 10 10 and 10 12 GeV are plotted. These show 
clear evidence that fluctuations become less important with increasing en- 
ergy. In figure 13 the maximum acoustic pulse height per GeV at z =8 m 
and r=l km from the shower is presented. For both the original and SVD 
parameterised data one hundred showers are modelled in each half decade. 
The output directly from Corsika is used to create the "direct" plot. The 
pulse height from interpolated SVD parameterisation are plotted at 1/4 and 
3/4 decades for clarity. The error-bars are drawn at the 10 and 90 centiles. 
The earlier ACORNE parameterisation [13] is shown for comparison. Neither 
of the parameterisations is perfect, however the differences are small; only a 
few percent at all energies. The SVD parameterisation model seems to slightly 
under-estimate the fluctuations. This is due to the statistics being slightly non 
Gaussian (leptokurtic). The earlier ACORNE parameterisation also works well 
and produces pulses with similar accuracies. 



5.5 Summary of the SVD parameterisation 



The SVD technique successfully models the shape of both the radial and lon- 
gitudinal distributions of the showers to a high accuracy (> 95% of the infor- 
mation is retained). However for the acoustic detection of cosmic ray neutrinos 
it is primarily energies of greater than 10 9 GeV that are of significance to the 
acoustic technique. In this region fluctuations only contribute a few percent to 
the acoustic pulse amplitude and our earlier parameterisation [15] is perfectly 
adequate for modelling pulses in this region. In broader terms however the 
SVD method outlined is a linear technique and does not rely on optimisation 
and as the accuracy is dictated by the number of eigenvectors, a trade off can 
be made between accuracy and complexity. 



6 Conclusions 



A new way of computing the acoustic signal for neutrino showers in water and 
ice has been described. The method is computationally fast and allows the 
most up to date knowledge of the attenuation to be incorporated naturally. 
This is now known to be complex in nature. A parameterisation of this at- 
tenuation is given. The properties of the expected acoustic signals from such 
showers have been described. These properties will be used in a future search 
for such interactions in an acoustic array. A matrix method of parameterising 
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the signals which includes fluctuations has also been described. This method 
is based on the SVD technique and it is shown to model both the shape of the 
radial and longitudinal distributions of the showers to a good accuracy. 
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b) Comparison of Triaxial and Spherical deposition 
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Fig. 1. Acoustic pulses from a 1 J deposition with an observer at 10 km, a) compar- 
ison of analytical and MC pulse for a spherical deposition a =1 cm b) comparison 
of a spherical [o =1 cm) and triaxial deposition [a x =1 cm, a y =0.5 cm, a z =1 m) 
c) rotating the observer from x-axis to the y-axis in the x — y plane for the triaxial 
distribution d) rotating the observer in the x — z plane from the x-axis to the z-axis 
for the triaxial distribution. 
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Sound Attenuation Coefficients 
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Fig. 2. Sound attenuations in sea water for the SAUND/Learned (/ 2 , Niess and 
Bertin, Francois and Garrison, Ainslie and McColm and the ACORNE parameter- 
isations. Since the Ainslie & McColm/ ACORNE and Francois & Garrison results 
match very closely they are depicted as the same curve in the figure. The attenu- 
ation in Antarctic ice is shown for comparison. Inset: the phase shifts at 1 km for 
the two complex attenuation models. 
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Pulse from a 6(x)6(t) deposition of 10 11 GeV of thermal energy at 1km 
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Fig. 3. The acoustic pulse for a point deposition of energy and an observer at 
lkm from the source for the f 2 , Niess and Bertin, Francois and Garrison, Ainslie 
and McColm and ACORNE parameterisations. The anticipated acoustic pulse for 
Antarctic ice is shown for comparison. The ACORNE and Niess and Bertin pulses 
are delayed due to the use of complex attenuation. 
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Fig. 4. The pulse at 1km for 10 5 — 10 12 GeV showers in water, a) pulse shape b) 
maximum and minimum pulse heights vs. energy 
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a) Sqrt of Pulse Energy 
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d) Captured Energy through a 1km Integrating Sphere 

5r 



> 
CD 



-10 



-15 




6 8 10 

Log 10 (E ) (GeV) 



12 



Fig. 5. Variation of acoustic pulse energy with angle in sea water a) square root of 
the pulse energy versus angle b) full width for 10% and 50% of maximum versus 
energy c) fluence through an integrating sphere d) captured energy versus deposited 
energy and linear fit (see equation 22). 
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Fig. 6. Mean frequency, Pressure x distance and asymmetry for the acoustic pulse 
in the pancake plane for a 10 11 GeV shower in sea water plotted as a function of 
distance. 
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c) Asymmetry d) Mean Frequency (kHz) 
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Fig. 7. Pulse parameters in sea water for a 10 11 GeV primary as a function of distance 
and angle a) maximum pulse height b) minimum pulse height c) pulse asymmetry 
d) mean frequency 
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a) Pulses at 1km from Cascade 
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Fig. 8. The pulse at 1km for 10 5 — 10 12 GeV showers in ice. a) pulse shape, b) 
maximum and minimum pulse heights vs. energy. 
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a) Sqrt of Pulse Energy b) Full Width of Pancake 
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Fig. 9. Variation of acoustic pulse energy with angle in ice a) square root of the 
pulse energy versus angle b) full width for 10% and 50% of maximum versus energy 
c) fluence through an integrating sphere d) captured energy versus deposited energy 
and linear fit (see equation 22). 
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Fig. 10. Pulse parameters in ice for a 10 GeV primary as a function of distance 
and angle a) maximum pulse height b) minimum pulse height c) pulse asymmetry 
d) mean frequency. 
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Fig. 11. Longitudinal distribution of four CORSIKA generated 10 5 GeV showers 
(points) with smoothing (dotted curve) and SVD reconstruction (solid curve). 
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Fig. 12. Fluctuation of acoustic pulses at 
and 10 12 GeV. 
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Comparison of Raw and Parameterised pulse hights (r=1km, z=8m) 
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Fig. 13. Pulse heights for CORSIKA generated showers at 1km from the source as a 
function of energy, with the SVD parameterisation and functional parameterisation 
shown for comparison. The 10th and 90th centiles are shown for the CORSIKA and 
SVD cases. 
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